﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using DMC=DiffusionMonteCarlo;

namespace Examples
{
    public class HydrogenAtom
    {
        private string result;
        public string Result
        {
            get { return result; }
        }
        public void Compute()
        {
            DMC.Input input = new DMC.Input();

            input.NumberOfParticles = 1;
            input.SpaceDimension = 3;

            input.Seed = 22312345L;
            input.InitialNumberOfReplicas = 500;
            input.RunTime = 1000.0;
            input.ValueOfTimeStep = 0.05;

            double[,] minmax = { {-10.0, 10.0},
                                 {-10.0, 10.0},
                                 {-10.0, 10.0} };
            input.LimitsOfCoordinates =  minmax;
            double[,] xyz = { {0.1, 0.1, 0.1} };
            input.InitialPosition = xyz;
            double[] delx = { 1.0, 1.0, 1.0 };
            input.ValueOfLengthInterval = delx;

            DMC.PotentialEnergy<double[,], double> V
                = new DMC.PotentialEnergy<double[,], double>(Function);
    
            DMC.PathIntegral pi 
                = new DMC.PathIntegral(input, V);
            
            DMC.Output output = pi.Compute();

            result += output.GroundStateEnergy.ToString() + "\r\n";
            for(int i = 0; i < output.GroundStateWaveFunction.Length; i++)
                result += output.GroundStateWaveFunction[i].ToString() + "\t";


            int[] nn = new int[input.SpaceDimension];
            for (int i = 0; i < input.SpaceDimension; i++)
            {
                nn[i] = (int)Math.Floor((input.LimitsOfCoordinates[i, 1] - input.LimitsOfCoordinates[i, 0])
                    / input.ValueOfLengthInterval[i]);
            }

            double[,,] Psi = new double[nn[0],nn[1],nn[2]]; 
            for (int i = 0; i < output.GroundStateWaveFunction.Length; i++)
            {
                int[] ii = MultiIndex(input.NumberOfParticles, nn, i);
                Psi[ii[0],ii[1],ii[2]] = output.GroundStateWaveFunction[i];
            }

            double[,] xy = new double[3, Psi.GetLength(0)];
            for (int i = 0; i < Psi.GetLength(0); i++)
            {
                xy[0, i] = input.LimitsOfCoordinates[0, 0]
                    + ((double)i+0.5) * input.ValueOfLengthInterval[0];
                xy[1, i] = output.GroundStateEnergy;
                xy[2, i] = Psi[i, nn[1]/2, nn[2]/2];
            }

            Graph graph = new Graph(xy);
            graph.Show();
        }
        
        private int[] MultiIndex(int np, int[] nn, int i)
        {
            int[] t = new int[np*nn.Length];
            int wh = np * nn.Length - 1;
          
            for (int j = np; j >= 1; j--)
            {
                for (int k = nn.Length - 1; k >= 0; k--)
                {
                    t[wh] = i % nn[k];
                    i -= t[wh];
                    i /= nn[k];
                    wh--;
                }
            }

            return t;
        }


        private double Function(double[,] x)
        {
            return -1.0 / Math.Sqrt(x[0, 0] * x[0, 0] + x[0, 1] * x[0, 1] + x[0, 2] * x[0, 2]);
        }
    }
}